SUBROUTINE READ_ECOCLIMAP &
& (dx,dy,xlat,xlon,imonth,nx,ny,coords_syst,xlake, &
& watdepth,alb,z0,fcveg,veg,fnmap)

! Subroutine READ_ECOCLIMAP reads ECOCLIMAP data
! on surface characteristics distribution


use MPI_VARIABLES, only : &
& nxs2i,nxs2e,nys2i,nys2e, &
& nx0i,nx1e,ny0i,ny1e, &
& rank_comm3d,comm3d,parallel
use ALLOC, only : &
& rearth

implicit none

include 'mpif.h'

include 'soil.inc'
include 'params_ECOCLIMAP.inc'

! Input/output variables

real(8), intent(inout) :: dx, dy 
real(8), intent(in) :: xlat, xlon ! Coordinates of domain center
integer(4), intent(in) :: imonth ! Currently one month is used, in future
                                 ! updating of monthly vegetation distribution
                                 ! during model integration should be done
integer(4), intent(in) :: nx, ny
integer(4), intent(in) :: coords_syst

real(8), intent(out) :: xlake(nxs2i:nxs2e,nys2i:nys2e)
real(8), intent(inout) :: watdepth(nxs2i:nxs2e,nys2i:nys2e)
real(8), intent(out) :: alb(nxs2i:nxs2e,nys2i:nys2e)
real(8), intent(out) :: z0(nxs2i:nxs2e,nys2i:nys2e)
real(8), intent(out) :: fcveg(1:nv2,nxs2i:nxs2e,nys2i:nys2e)
real(8), intent(out) :: veg(nxs2i:nxs2e,nys2i:nys2e)

character(len=80), intent(in) :: fnmap

! Local variables

real(8), parameter :: pi = 3.141592
real(8), parameter :: z0_water = 1.d-3, alb_water = 5.d-2
real(8), parameter :: watdepthref = 5. ! Reference lakes and seas depth, m
real(8) :: zmaxlon_m, zminlon_m, zmaxlat_m, zminlat_m
real(8) :: dlat_ECO, dlon_ECO, dx_ECO, dy_ECO, x_ECO, x_m
real(8) :: halfsize
real(8) :: xx(1:2), xxx
real(8), allocatable :: workdp(:,:), work2(:,:), workdp2(:,:)
integer(4), parameter :: iunit = 1, irecl = 1 ! The length of record for one real value, 
                                              ! for Intel compiler is measured in 4-byte words!
real, parameter :: codeabs = 999. ! In ECOCLIMAP - code of data absence, 
                                  ! indicates sea points
real, allocatable :: work(:,:), work1d(:)
integer(4) :: res_km, nrec, nx_ECO, ny_ECO, ixstart, iystart, n_ECO
integer(4) :: i, j, k, m
integer(4) :: ierr
character(len=50) :: filename
character(len=80) :: filename_out
character(len=20) :: filenames(1:nv2)
character :: ext_res*3, ext_month*2, workch*2

character(len=80), external :: FEXT

!--------------------------------------------------------------------------
!VEGEATION TYPES for STOMATA from ECOCLIMAP ARE
!1-broadleaf evergreen tree
!2-broadleaf decidious tree
!3-needleleaf tree
!4-grass
!5-shrub
!6-arctic shrub
!7-crop
!--------------------------------------------------------------------------

if (coords_syst == 1) then

  res_km = nint(dx/1000.)
  if (res_km /= nint(dy/1000.)) then
    call OPTWARNSTOP('In READ_ECOCLIMAP: dx /= dy!: STOP')
  endif

  halfsize = 0.5*(nx-1)*dx/(rearth*dcos(xlat*pi/180.))/pi*180. ! degrees
  zmaxlon_m = xlon + halfsize
  zminlon_m = xlon - halfsize

  halfsize = 0.5*(ny-1)*dy/rearth/pi*180. ! degrees
  zmaxlat_m = xlat + halfsize
  zminlat_m = xlat - halfsize

elseif (coords_syst == 2) then
  zminlon_m = zminlon_ECO
  zmaxlon_m = zmaxlon_ECO
  zminlat_m = zminlat_ECO
  zmaxlat_m = zmaxlat_ECO
endif

!if (rank_comm3d == 0) write(*,*) ' zminlon_m = ', zminlon_m, &
!&                                ' zmaxlon_m = ', zmaxlon_m, &
!&                                ' zminlat_m = ', zminlat_m, &
!&                                ' zmaxlat_m = ', zmaxlat_m, &
!&                                ' zminlon_ECO = ', zminlon_ECO, &
!&                                ' zmaxlon_ECO = ', zmaxlon_ECO, &
!&                                ' zminlat_ECO = ', zminlat_ECO, &
!&                                ' zmaxlat_ECO = ', zmaxlat_ECO


if (.not. &
  & (zminlon_m >= zminlon_ECO .and. zmaxlon_m <= zmaxlon_ECO .and. &
  &  zminlat_m >= zminlat_ECO .and. zmaxlat_m <= zmaxlat_ECO) ) then
  call OPTWARNSTOP('In READ_ECOCLIMAP: model domain exceeds ECOCLIMAP data &
  & domain')
endif

nx_ECO = nint((zmaxlon_ECO-zminlon_ECO)*float(nres_ECO)/float(nscale_ECO))
ny_ECO = nint((zmaxlat_ECO-zminlat_ECO)*float(nres_ECO)/float(nscale_ECO))
nrec = nx_ECO ! Record length

if (coords_syst == 2) then
  if (nx /= nx_ECO .or. ny /= ny_ECO) then
    call OPTWARNSTOP('In READ_ECOCLIMAP: nx /= nx_ECO or ny /= ny_ECO while &
    & spherical coordinates are used')
  endif
endif

! nres_ECO - number of grid steps of ECOCLIMAP in 1 degree latitude and longitude
dlon_ECO = 1./float(nres_ECO)*nscale_ECO
dlat_ECO = 1./float(nres_ECO)*nscale_ECO
dx_ECO = dlon_ECO*(rearth*dcos(xlat*pi/180.))*pi/180.
dy_ECO = dlat_ECO*pi/180.*rearth

if (coords_syst == 2) then
! If spherical coordinates are used, dx and dy are determined by
! latitude and longitude grid increments and latitude of the center 
! of domain
  dx = dx_ECO
  dy = dy_ECO
endif

iystart = nint((zmaxlat_ECO - zmaxlat_m)/dlat_ECO) + 1
!iystart = nint((zminlat_m - zminlat_ECO)/dlat_ECO)
ixstart = nint((zminlon_m - zminlon_ECO)/dlon_ECO) + 1

!if (rank_comm3d == 0) write(*,*) ' nx_ECO = ', nx_ECO, &
!&                                ' ny_ECO = ', ny_ECO, &
!&                                ' ixstart = ', ixstart, &
!&                                ' iystart = ', iystart

!write(ext_res,'(i3.3)') res_km
write(ext_res,'(i3.3)') nscale_ECO
write(ext_month,'(i2.2)') imonth

filenames(1) = 'f_tropical_tree'
filenames(2) = 'f_broadleaf_tree'
filenames(3) = 'f_conif_tree'
filenames(4) = 'f_grassland'
filenames(5) = 'f_park_marshes'
filenames(6) = 'f_trop_grass'
filenames(7) = 'f_C3_crops'

! Reading vegetation types coverage
allocate (work(0:nx+1,0:ny+1))
allocate (work2(0:nx+1,0:ny+1), &
& workdp(0:nx+1,0:ny+1), workdp2(0:nx+1,0:ny+1))
work2 = 0. ! Coverage by all vegetation types
do k = 1, nv2
  filename = 'map/ECOCLIMAP/'//filenames(k)(1:len_trim(filenames(k)))// &
  & '.'//ext_res
  if (rank_comm3d == 0) then
    call READFILE
    call FILLABS(0.d0)
  endif
  if (parallel) then
    call MPI_BCAST(work,(nx+2)*(ny+2),MPI_REAL,0,comm3d,ierr)
  endif
  workdp = work
  if (rank_comm3d == 0) then
    write (workch,'(i2.2)') k
    filename_out = FEXT(fnmap,'v'//workch)
    call WRIGRD(workdp, &
    &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
    &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out)
  endif
  call EXTRAH_MPI &
  & (workdp, 0, nx+1, 0, ny+1, 0, 0, &
  & 0, nx+1, 0, ny+1, 0, 0)
  fcveg(k,nx0i:nx1e,ny0i:ny1e) = workdp(nx0i:nx1e,ny0i:ny1e)
  work2 = work2 + workdp
enddo
veg(nx0i:nx1e,ny0i:ny1e) = work2(nx0i:nx1e,ny0i:ny1e)
if (rank_comm3d == 0) then
  filename_out = FEXT(fnmap,'veg')
  call WRIGRD(work2, &
  &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
  &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out)
endif

! Reading water bodies coverage
filename = 'map/ECOCLIMAP/F_wat.'//ext_res
if (rank_comm3d == 0) then
  call READFILE
  call FILLABS(0.d0)
endif
if (parallel) then
  call MPI_BCAST(work,(nx+2)*(ny+2),MPI_REAL,0,comm3d,ierr)
endif
workdp = work
do j = 0, ny+1
  do i = 0, nx+1
    if (workdp(i,j) > 0.d0) then
      workdp2(i,j) = watdepthref
    endif
  enddo
enddo
if (rank_comm3d == 0) then
  filename_out = FEXT(fnmap,'lak')
  call WRIGRD(workdp, &
  &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
  &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out)
  filename_out = FEXT(fnmap,'dep')
  call WRIGRD(workdp2, &
  &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
  &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out)
endif
call EXTRAH_MPI &
& (workdp, 0, nx+1, 0, ny+1, 0, 0, &
& 0, nx+1, 0, ny+1, 0, 0)
xlake(nx0i:nx1e,ny0i:ny1e) = workdp(nx0i:nx1e,ny0i:ny1e)
call EXTRAH_MPI &
& (workdp2, 0, nx+1, 0, ny+1, 0, 0, &
& 0, nx+1, 0, ny+1, 0, 0)
watdepth(nx0i:nx1e,ny0i:ny1e) = workdp2(nx0i:nx1e,ny0i:ny1e)

! Reading albedo
filename = 'map/ECOCLIMAP/alb.'//ext_month//'.'//ext_res
if (rank_comm3d == 0) then
  call READFILE
  call FILLABS(alb_water)
endif
if (parallel) then
  call MPI_BCAST(work,(nx+2)*(ny+2),MPI_REAL,0,comm3d,ierr)
endif
workdp = work
if (rank_comm3d == 0) then
  filename_out = FEXT(fnmap,'alb')
  call WRIGRD(workdp, &
  &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
  &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out) 
endif
call EXTRAH_MPI &
& (workdp, 0, nx+1, 0, ny+1, 0, 0, &
& 0, nx+1, 0, ny+1, 0, 0)
alb(nx0i:nx1e,ny0i:ny1e) = workdp(nx0i:nx1e,ny0i:ny1e)


! Reading roughness
filename = 'map/ECOCLIMAP/z0.'//ext_month//'.'//ext_res
if (rank_comm3d == 0) then
  call READFILE
  call FILLABS(z0_water)
endif
if (parallel) then
   call MPI_BCAST(work,(nx+2)*(ny+2),MPI_REAL,0,comm3d,ierr)
endif
workdp = work
if (rank_comm3d == 0) then
  filename_out = FEXT(fnmap,'z0 ')
  call WRIGRD(workdp, &
  &  0.d0, (nx-1)*dx/1000., 0.d0, (ny-1)*dy/1000., xx(1), xx(2), &
  &  0,nx+1,0,ny+1,0,nx+1,0,ny+1,filename_out)
endif
call EXTRAH_MPI &
& (workdp, 0, nx+1, 0, ny+1, 0, 0, &
& 0, nx+1, 0, ny+1, 0, 0)
z0(nx0i:nx1e,ny0i:ny1e) = workdp(nx0i:nx1e,ny0i:ny1e)

! Check the coverages for consistency: bare soil coverage >= 0
!do j = ny0i, ny1e
!  do i = nx0i, nx1e
!    if (xlake(i,j) > 1.) then
!      write(*,*) 'veg(i,j) = ', veg(i,j), 'xlake(i,j) = ', xlake(i,j)
!      call MPI_BARRIER(comm3d,ierr)
!      call OPTWARNSTOP('In READ_ECOCLIMAP: inconsistent vegetation &
!      & and water bodies coverage: STOP')
!    endif
!    do k = 1, nv2
!      if (fcveg(k,i,j) == codeabs) fcveg(k,i,j) = 0.
!    enddo
!    if (veg(i,j) == nv2*codeabs) veg(i,j) = 0.
!    if (z0(i,j) == codeabs) z0(i,j) = z0_water
!    if (alb(i,j) == codeabs) alb(i,j) = alb_water
!  enddo
!enddo

deallocate (work)
deallocate (work2,workdp,workdp2)
contains


SUBROUTINE READFILE

!INQUIRE(file=filename, recl=nrec)
open(unit=iunit,file=filename,access='direct', &
& form='unformatted',recl=nrec*irecl,status='old')
work = 0.
allocate(work1d(1:nx_ECO))

!if (filename(1:len_trim(filename)) == 'map/F_wat.010') then
!  read(unit=iunit,rec=60) work1d(:)
!  do j = 1, nx_ECO
!    write(*,*) work1d(j)
!  enddo
!endif

do j = iystart, iystart + ny - 1
  read(unit=iunit,rec=j) work1d(:)
! Interpolation is implemented only for x coordinate,
! in y coordinate spacing of ECOCLIMAP data is assumed
! to be equal to that of the model
  if (coords_syst == 1) then
    work(1,iystart+ny-j) = work1d(ixstart)
    x_ECO = 0.
    x_m = 0.
    n_ECO = ixstart
    do m = 2, nx
      x_m = x_m + dx
      do while (x_ECO < x_m)
        x_ECO = x_ECO + dx_ECO
        n_ECO = n_ECO + 1
      enddo
      xxx = (x_ECO-x_m)/dx_ECO
      if (work1d(n_ECO-1) == codeabs .or. &
      &   work1d(n_ECO) == codeabs) then
        work(m,iystart+ny-j) = codeabs
      else
        work(m,iystart+ny-j) = work1d(n_ECO-1)*xxx + work1d(n_ECO)*(1.-xxx)
      endif
    enddo
  elseif (coords_syst == 2) then
    work(1:nx,iystart+ny-j) = work1d(ixstart:ixstart+nx-1)
  endif

enddo
deallocate(work1d)
close(iunit)

END SUBROUTINE READFILE


SUBROUTINE FILLABS(fillvalue)

real(8), intent(in) :: fillvalue

do j = 1, ny
  do i = 1, nx
    if (work(i,j) == codeabs) work(i,j) = fillvalue
  enddo
enddo


END SUBROUTINE FILLABS


END SUBROUTINE READ_ECOCLIMAP
